Method for processing a digital image

ABSTRACT

A method and system for processing a distorted digital image B that is a convolution of an undistorted image F and a point spread function. Noise is removed from the image B so as to produce an image B′ of reduced noise. The image F is then obtained based upon a calculation involving the image B′.

FIELD OF THE INVENTION

This invention relates to methods for processing a digital image.

BACKGROUND OF THE INVENTION

Global distortion, or blurring, of a picture can arise from variousfactors such as distortion due to out-of-focus optics and lineartranslation or shaking of the camera during the exposure time.

Blurring of a digital image may be described by means of a convolution:B ₀(x)=∫dx′F(x′)h(x−x′)  (1)

where B₀(x) is the intensity of the pixel at the address x=(x,y) in thedistorted picture, x being a two-dimensional vector, F(x) is theintensity of the pixel x in the undistorted image, and h(x) is theso-called point spread function (PSF) that describes the distortion. Thefunction B₀(x) is typically obtained as the output from a digitalcamera. The PSF for an image distorted by out-of-focus optics, forexample, can be described in a first approximation by a function h thatis constant inside a circle of radius r and h(x)=0 outside the circle.

Rectifying a distorted image involves determining the function F giventhe functions B₀ and h. The convolution (1) can be Fourier transformedto yield{tilde over (B)} ₀(q)={tilde over (F)}(q)·{tilde over (h)}(q)  (2)where “{tilde over ( )}” represents the Fourier transform. Hence,

$\begin{matrix}{{\overset{\sim}{F}(q)} = \frac{{\overset{\sim}{B}}_{0}(q)}{\overset{\sim}{h}(q)}} & (3)\end{matrix}$In principle, therefore, F(x) may be obtained from the inverse Fouriertransform of {tilde over (B)}(q)/{tilde over (h)}(q). In practice,however, this solution is characterized by a very low signal-to-noiseratio (SNR), due to amplification of noise at frequencies q at which{tilde over (h)}(q) is very small.

SUMMARY OF THE INVENTION

In the following description and set of claims, two explicitlydescribed, calculable or measurable variables are considered equivalentto each other when the two variables are substantially proportional toeach other.

The present invention provides a method for rectifying a distorteddigital image B₀(x) to produce a rectified image F(x). In accordancewith the invention, noise is removed from the function {tilde over(B)}₀(q) before applying the inverse Fourier transform to the right sideof equation. (3). A noise function N is used to evaluate the amount ofnoise for functions {tilde over (B)}(q) that deviate slightly from{tilde over (B)}₀(q) and a new function {tilde over (B)}(q) is selectedthat deviates slightly from {tilde over (B)}₀(q) and that minimizes thenoise function N. In a preferred embodiment of the invention, the noiseN in an image B(x) is obtained based upon a calculation involving thegradient of the function P obtained by inverse Fourier transform of{tilde over (B)}(q)/{tilde over (h)}(q). In a most preferred embodiment,the noise N is calculated according to the equation:N=∫∇P(x)·∇P′(x)dx  (4)where “*” indicates complex conjugate. Equation (4) may be written inthe equivalent formN=∫dq·q ² ∥{tilde over (D)}(q)∥² ·∥{tilde over (B)} ₀(q)∥²  (4)where D(x) is the so-called deconvolution filter (DCF) defined as1/h(x), wherein h(x) is the point spread function characteristic of thedistortion, and q is a two-dimensional vector in Fourier space. Therectified image F(x) may then be obtained by calculating the inverseFourier transform of the right side of equation (3) using the thusobtained B(x). However, a pattern characteristic of the DCF may havebeen superimposed on the obtained F(x). This pattern originates in theDCF and is correlated with it. In the most preferred embodiment, thispattern is removed.

The invention thus concerns a method for processing a digital image B₁,the image B₁ being a convolution of an image F and a point spreadfunction h, comprising removing noise from the image B₁ so as to producean image B′ of reduced noise, and calculating F based upon B′.

The invention further concerns a method for processing a deconvolutedimage B, the image B having been deconvoluted according to adeconvolution filter D, the method comprising reducing correlationbetween the image and the deconvolution filter.

Yet still further the invention concerns a method for processing adigital image B₁, the image B₁ being a convolution of an image F and apoint spread function h comprising the steps of:

-   -   (a) removing noise from the image B₁ so as to produce an image        B′ of reduced noise;    -   (b) obtaining function {tilde over (P)}₁(q) according to the        algebraic expression {tilde over (P)}₁(q)={tilde over        (B)}′(q)/{tilde over (h)}(q);    -   (c) reducing correlation between {tilde over (P)}₁ and 1/{tilde        over (h)} so as to product a function {tilde over (P)}′ of        reduced correlation; and    -   obtaining a rectified image F by inverse Fourier transform of        {tilde over (P)}′(q).

By yet still another aspect the invention concerns a method forobtaining a radius r of a point spread function h describing anout-of-focus distortion of a digital image B, the method comprising astep of calculating a gradient at a plurality of pixels in the image B.

By a further aspect the invention concerns a program storage devicereadable by machine, tangibly embodying a program of instructionsexecutable by the machine to perform method steps for processing adigital image B₁, the image B₁ being a convolution of an image F and apoint spread function h, comprising removing noise from the image B₁ soas to produce an image B′ of reduced noise, and calculating F based uponB′.

Further the invention concerns a program storage device readable bymachine, tangibly embodying a program of instructions executable by themachine to perform method steps for processing a deconvoluted image B,the image B having been deconvoluted according to a deconvolution filterD, the method comprising reducing correlation between the image and thedeconvolution filter.

Yet still further the invention concerns a program storage devicereadable by machine, tangibly embodying a program of instructionsexecutable by the machine to perform method steps for processing adigital image B₁, the image B₁ being a convolution of an image F and apoint spread function h comprising the steps of.

-   -   (a) removing noise from the image B₁ so as to produce an image        B′ of reduced noise;    -   (b) obtaining function {tilde over (P)}₁(q) according to the        algebraic expression {tilde over (P)}₁(q)={tilde over        (B)}′(q)/{tilde over (h)}(q);    -   (c) reducing correlation between {tilde over (P)}₁ and {tilde        over (1)}/h so as to product a function {tilde over (P)}′ of        reduced correlation; and    -   (d) obtaining a rectified image F by inverse Fourier transform        of {tilde over (P)}′(q).

Further the invention concerns a program storage device readable bymachine, tangibly embodying a program of instructions executable by themachine to perform method steps for obtaining a radius r of a pointspread function h describing an out-of-focus distortion of a digitalimage B, the method comprising a step of calculating a gradient at aplurality of pixels in the image B.

BRIEF DESCRIPTION OF THE DRAWINGS

In order to understand the invention and to see how it maybe carried outin practice, a preferred embodiment will now be described, by way ofnon-limiting example only, with reference to the accompanying drawings,in which:

FIG. 1 shows a digital image showing distortion due to out of focusoptics;

FIG. 2 shows a histogram of heights of steps in the image of FIG. 1;

FIG. 3 shows the result of rectifying the image of FIG. 1 in accordancewith the invention; and

FIG. 4 shows a system for digital image processing in accordance with anembodiment of the invention.

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT

In a preferred embodiment, as shown in FIG. 4, a system 20 for digitalimage processing comprises a camera 22 and an image processor 24 (whichmay be a separate unit as shown or may be integrated in the camera). Thecamera typically comprises objective optics 26, having a certain pointspread function (PSF). Optics 26 focus an image onto an image sensor 28,comprising multiple pixels 30. Image sensor 28 outputs an image that isdistorted, for reasons explained below. Digital image processor 24processes the image to remove noise and distortion. The image processormay comprise a programmable processor, which operates in accordance withprogram instructions stored in a program storage device 32. Prior toremoving noise from the distorted image B₀(x), the function B₀(x) ispreprocessed as follows. Firstly, if the output signal B(x) from thedigital camera is not linear with the input intensities I(x), the imageB₀(x) is transformed to make the signal linear with the intensity. Thisis accomplished by applying to the signal B₀(x) the inversetransformation that was applied by the camera to the intensities inorder to produce the signal B(x). For example, if the camera performs atransformation on the intensity of the form B(x)=A(A⁻¹I(x))^(γ), where Ais a scaling factor and γ a fixed exponent, the signal B₀(x) istransformed by A(B₀(x)/A)^(γ) ⁻¹ to obtain a new B(x). For many digitalcameras, γ=⅓ in order to make the obtained image more pleasing to theeye. The function B(x) is then transformed so as to decay smoothly atthe edges to zero in order to make the image periodic at the edges. Thefunction B₁(x) produced by this preprocessing is operated on inaccordance with the invention as described in detail below.

In accordance with the invention, a noise function N is used to evaluatethe amount of noise for functions {tilde over (B)}(q) that deviateslightly from {tilde over (B)}₁(q) and a new function {tilde over(B)}(q) is selected that deviates slightly from {tilde over (B)}₁(q) andthat minimizes the noise function N. In a preferred embodiment of theinvention, the noise N in an image {tilde over (B)}(q) is obtained basedupon a calculation involving the gradient of the function P obtained byinverse Fourier transform of {tilde over (B)}(q)/{tilde over (h)}(q). Ina most preferred embodiment, the noise N is calculated according to theequation:N=∫∇P(x)·∇*(x)dxor equivalently, N=∫dq·q ² ∥{tilde over (D)}(q)∥² ·∥{tilde over (B)}(q)∥²  (4)

where D(x) is the deconvolution filter defined as 1/h (x), where h(x) isthe point spread function characteristic of the distortion.

A function B(x) of essentially minimal noise that deviates only slightlyfrom B₁(x) may be found by evaluating the variation of N with respect toB*(q), (∂N/∂{tilde over (B)}*(q)). For example, a sequence of functions{tilde over (B)}₁(q) may be generated iteratively by:

${{\overset{\sim}{B}}_{l + 1}(q)} = {{{\overset{\sim}{B}}_{l}(q)} - {ɛ{\frac{\partial N}{{\partial{\overset{\sim}{B}}_{i}}*(q)}.}}}$

In the most preferred embodiment, N is given by (4) so that

$\begin{matrix}{{{\overset{\sim}{B}}_{l + 1}(q)} = {{{{\overset{\sim}{B}}_{i}(q)} - {ɛ\frac{\partial N}{{\partial{\overset{\sim}{B}}_{l}}*(q)}}} = {{{{\overset{\sim}{B}}_{l}(q)} - {ɛ{{{\overset{\sim}{B}}_{i}(q)} \cdot {{\overset{\sim}{D}(q)}}^{2} \cdot q^{2}}}} = {{{\overset{\sim}{B}}_{i}(q)}{\left( {1 - {ɛ{{\overset{\sim}{D}(q)}}^{2}q^{2}}} \right) \cdot}}}}} & (5)\end{matrix}$

Thus in the most preferred embodiment, (N given by (4)), for finite α,i=n for large n, and ε=α/n, {tilde over (B)}_(n)(q) is given by:

${{{\overset{\sim}{B}}_{n}(q)} = {{{\overset{\sim}{B}}_{1}(q)}\left( {1 - {\frac{\alpha}{n}{{\overset{\sim}{D}(q)}}^{2}q^{2}}} \right)^{n}}},$

and in the limit,

$\begin{matrix}{{{\overset{\sim}{B}}_{\alpha}(q)} = {{\lim\limits_{n\rightarrow\infty}{{{\overset{\sim}{B}}_{1}(q)}\left( {1 - {\frac{\alpha}{n}{{\overset{\sim}{D}(q)}}^{2}q^{2}}} \right)^{n}}} = {{{\overset{\sim}{B}}_{1}(q)}{\mathbb{e}}^{{- \alpha}{{\overset{\sim}{D}{(q)}}}^{2}q^{2}}}}} & (6)\end{matrix}$

{tilde over (B)}_(α)(q) is preferably multiplied by a factor so that itacquires the same norm as {tilde over (B)}₁(q) to produce a function{tilde over (B)}(q) and a new function {tilde over (P)}₁(q)={tilde over(B)}(q)/{tilde over (h)}(q) is then obtained.

F(x) may now be obtained by inverse Fourier transform of {tilde over(P)}₁(q). However, a pattern characteristic of the DCF may have beensuperimposed on {tilde over (P)}₁(q). This pattern originates in the DCFand is correlated with it. In the most preferred embodiment, thispattern is removed from the function {tilde over (P)}₁(q) beforeapplying the inverse Fourier transform and return to x-space.

In order to remove this pattern, the correlation between the pattern andthe DCF is decreased, An overall correlation function C is used toevaluate the correlation between the DCP and {tilde over (P)}(q) forfunctions {tilde over (P)}(q) that deviate slightly from {tilde over(P)}₁(q) and a new function {tilde over (P)}(q) is selected thatdeviates slightly from {tilde over (P)}₁(q) and that minimizes thecorrelation function C. In a preferred embodiment of the invention, thecorrelation C between the DCF and a function {tilde over (P)}(q) iscalculated according to the equation:C=∫dq∥{tilde over (D)}(q)∥² ·∥{tilde over (P)}(q)∥²  (7)

The variation of C with respect to {tilde over (P)}*(q), ∂C/∂{tilde over(P)}* is used to change {tilde over (P)}(q) in order to reduce C. forexample, a sequence of functions {tilde over (P)}i (q) may be generatediteratively by:

${{\overset{\sim}{P}i} + {1(q)}} = {{\overset{\sim}{P}{i(q)}} - {ɛ{\frac{\partial C}{\partial{{\overset{\sim}{P}}^{*}(q)}}.}}}$

In the most preferred embodiment, C is given by (7) so that

${{{\overset{\sim}{P}}_{i + 1}(q)} = {{{{\overset{\sim}{P}}_{i}(q)} - {ɛ\frac{\partial C}{\partial{{\overset{\sim}{P}}^{*}(q)}}}} = {{{{\overset{\sim}{P}}_{i}(q)} - {ɛ \cdot {{\overset{\sim}{P}}_{i}(q)} \cdot {{\overset{\sim}{D}(q)}}^{2}}} = {{{\overset{\sim}{P}}_{i}(q)}\left( {1 - {ɛ{{\overset{\sim}{D}(q)}}^{2}}} \right)}}}},$

Thus, in the most preferred embodiment (C given by (7)), for finite β,i=n for large n, and ε=β/n,

$\begin{matrix}{{{\overset{\sim}{P}}_{n + 1}(q)} = {{{\overset{\sim}{P}}_{1}(q)}\left( {1 - {\frac{\beta}{n}{{\overset{\sim}{D}(q)}}^{2}}} \right)^{n}}} & (8)\end{matrix}$

and in the limit,

$\begin{matrix}{{{\overset{\sim}{P}}_{\beta}(q)} = {{\lim\limits_{n\rightarrow\infty}{{{\overset{\sim}{P}}_{1}(q)}\left( {1 - {\frac{\beta}{n}{{\overset{\sim}{D}(q)}}^{2}}} \right)^{n}}} = {{{\overset{\sim}{P}}_{1}(q)}{\mathbb{e}}^{{- \beta}{{\overset{\sim}{D}{(q)}}}^{2}}}}} & (9)\end{matrix}$

{tilde over (P)}_(β)(q) is multiplied by a factor to produce a {tildeover (P)}(q) having the same norm as {tilde over (P)}₁(q)·F(x) is thenobtained by inverse Fourier transfer of {tilde over (P)}(q).

After removing the superimposed pattern correction, the function F(x)may be converted from the linear intensity range to an intensity spacesuited for the eye. However, where the original B₀(x) was near 0, theslope of B₀(x) was very steep so that when values of B₀(x) near 0 areconverted for discrete representation, much information is lost. Thus,despite the fact that the resolution of the data is nominally 8 bit,with dark colors (corresponding to low values of B), the resolution ofthe data is equivalent to 3 bit.

Because much information has been lost in the dark regions, convertingF(x) to an intensity space suited for the eye is preferably notperformed simply by applying the inverse of the transformationpreviously applied by the camera to B(x) in order to convert it to thelinear range. Instead, image F(x) is preferably first mapped into aninterval, and then the inverse transformation is applied Thus, forexample, if the transformation A(B(x)/A)^(γ) ⁻¹ was applied to B(x) bythe camera, F(x) is mapped linearly into the interval [0,255] and thetransformation A(F(x)/A)^(γ) is applied to the function F(x). For theobtained image, a histogram is obtained from which a cut-off value isobtained at which a predetermined fraction of pixels, for example 10% ofthe pixels in the image have an intensity value less than or equal tothis cut-off value. The intensities are then linearly remapped so that255 is mapped to 255, the cut-off value is mapped to 25 (10% of 255),and intensity values between 0 and the cut-off are mapped to zero.

The PSF for an image distorted by out-of-focus optics can be describedin a first approximation by a function h that is constant inside acircle of radius r and h(x)=0 outside the circle. The radius r may bedetermined from a distorted image B(x) by the following algorithm. Thealgorithm makes use of the fact that in an unfocused picture, boundaries(referred to herein as “steps”) are gradual and not abrupt.

A step parallel to the y-axis at x=x₀ of an image is described by aHeaviside function Θ(x−x₀) independent of y. The mathematicaldescription of the removal from focus of this step is then given by:I(x)=∫dx′dy′h(x−x′,y−y′)Θ(x′−x ₀)where the intensity I in the vicinity of the step is independent of yand h(x,y) is the PSF function. Integration with respect to y yields

${I(x)} = {{\int{{\mathbb{d}x^{\prime}}{\Theta\left( {x^{\prime} - x_{0}} \right)}{\int{{\mathbb{d}\left( {y - y^{\prime}} \right)}{h\left( {{x - x^{\prime}},{y - y^{\prime}}} \right)}}}}} = {\int{{\mathbb{d}x^{\prime}}{\Theta\left( {x^{\prime} - x_{0}} \right)}{\frac{2\sqrt{r^{2} - \left( {x - x^{\prime}} \right)^{2}}}{\pi\; r^{2}}.}}}}$The slope of the step at x=x₀ is

$\begin{matrix}{{{\frac{{\mathbb{d}1}(x)}{\mathbb{d}x}}_{x = x_{0}} = {\frac{\mathbb{d}}{\mathbb{d}x}{\int{{\mathbb{d}x^{\prime}}{\Theta\left( {x^{\prime} - x_{0}} \right)}\frac{2\sqrt{r^{2} - \left( {x - x^{\prime}} \right)^{2}}}{\pi\; r^{2}}}}}}}_{x = x_{0}} \\{{{= {\int{{\mathbb{d}x^{\prime}}{\Theta\left( {x^{\prime} - x_{0}} \right)}\frac{\mathbb{d}}{\mathbb{d}x}\left( \frac{2\sqrt{r^{2} - \left( {x - x^{\prime}} \right)^{2}}}{\pi\; r^{2}} \right)}}}}_{x = x_{0}} =} \\{= {- {\int{{\mathbb{d}x^{\prime}}{\Theta\left( {x^{\prime} - x_{0}} \right)}\frac{\mathbb{d}}{\mathbb{d}x^{\prime}}\left( \frac{2\sqrt{r^{2} - \left( {x_{0} - x^{\prime}} \right)^{2}}}{\pi\; r^{2}} \right)}}}} \\{= {\int{{\mathbb{d}{x^{\prime}\left( \frac{\mathbb{d}{\Theta\left( {x^{\prime} - x_{0}} \right)}}{\mathbb{d}x^{\prime}} \right)}}\frac{2\sqrt{r^{2} - \left( {x_{0} - x^{\prime}} \right)^{2}}}{\pi\; r^{2}}}}} \\{= {{\int{{\mathbb{d}x^{\prime}}{\delta\left( {x^{\prime} - x_{0}} \right)}\frac{2\sqrt{r^{2} - \left( {x_{0} - x} \right)^{2}}}{\pi\; r^{2}}}} = {\frac{2r}{\pi\; r^{2}} = \frac{2}{\pi\; r}}}}\end{matrix}$In accordance with the invention, r is found by determining the slopes(x) at each pixel x where a step exists in the image B₁(x), where s isdefined as the absolute value of the gradient of B₁(x) divided by theheight of the step at x. The radius r(x) at x is

${r(x)} = \frac{2}{\pi\;{s(x)}}$

Edges in the image B₁(x) are detected by any method of edge detection asis known in the art, for example, as disclosed in Crane, R., “SimplifiedApproach to Image Processing, A: Classical and Modern Techniques in C”;Chapter 3. Prentice Hall, 1996, Calculated radii r(x) for pixels x at anedge in the image B₁(x) are normalized by dividing by the height of thestep at the edge. The normalized radii are arranged in a histogram. Theradius of the PSF is obtained essentially equal to the radius of maximumfrequency in the histogram.

EXAMPLE

FIG. 1 shows an image that is blurred due to out-of focus optics. Thisimage was rectified in accordance with the invention as follows.

The radius of the PSF of the distortion was first found as follows.Edges in the image were detected and radius r(x) was calculated for eachpixel x located at an edge as described above. Each calculated radiuswas normalized by dividing it by the height of the step at x. FIG. 2shows a histogram of the normalized heights. As can be seen in FIG. 2,the radius of maximum frequency of the histogram was found to be 4pixels, and this was taken as the radius of the PSF of the distortion.

{tilde over (B)}_(α)(q) was then obtained according to equation (6)using α=0.001, and then normalized so as to have the same norm as {tildeover (B)}₁(q), A function {tilde over (P)}₁(q) was then obtained byapplying the right side of equation (3) using {tilde over (B)}_(α)(q)for {tilde over (B)}₀(q)·{tilde over (P)}_(β)(q) was then obtainedaccording to equation (9) using β=0.01, and then normalized so as tohave the same norm as {tilde over (P)}₁(q). The function F(x) was thenobtained by inverse Fourier transform of {tilde over (P)}₁(q). Therectified image F(x) is shown in FIG. 3.

It will also be understood that the system according to the inventionmay be a suitably programmed computer. Likewise, the inventioncontemplates a computer program being readable by a computer forexecuting the method of the invention. The invention furthercontemplates a machine-readable memory tangibly embodying a program ofinstructions executable by the machine for executing the method of theinvention.

1. A method for processing a digital image B₁, the image B₁ being aconvolution of an image F and a point spread function h, comprisingremoving noise from the image B₁ so as to produce an image B′ of reducednoise, and calculating F based upon B′, wherein an amount of noise iscalculated in a plurality of images B, and the image B′ is selected asan image of essentially minimal noise among the images B.
 2. A methodfor processing a digital image B₁, the image B₁ being a convolution ofan image F and a point spread function h, comprising removing noise fromthe image B₁ so as to produce an image B′ of reduced noise, andcalculating F based upon B′, wherein the amount of noise in an image iscalculated using an algebraic expression involving the gradient of afunction P(x) obtained by inverse Fourier transform of {tilde over(B)}(q)/{tilde over (h)}(q).
 3. The method of claim 2, wherein theamount of noise N in an image B is calculated according to the algebraicexpression N=∫∇P(x)·∇P*(x)dx, wherein Δ indicates the gradient and “*”indicates complex conjugate.
 4. The method according to claim 3 wherein{tilde over (B)}′(q), the Fourier transform of B′, is equal to {tildeover (B)}_(i+1)(q) for some integer i, where {tilde over (B)}_(i+1)(q)is obtained according to the algebraic expression {tilde over(B)}_(i+1)(q)={tilde over (B)}₁(q)(1+ε∥{tilde over (D)}(q)∥² q ²)^(i),where ε is a small positive number.
 5. The method according to claim 3wherein {tilde over (B)}′(q) is obtained according to the algebraicexpression {tilde over (B)}′(q)={tilde over(B)}₁(q)e^(−α∥{tilde over (D)}(q)∥) ² ^(q) ² , where α is apredetermined constant, and {tilde over (D)}(q) is the Fourier transformof 1/h.
 6. A method for processing a digital image B₁, the image B₁being a convolution of an image F and a point spread function h,comprising removing noise from the image B₁ so as to produce an image B′of reduced noise, and calculating F based upon B′, wherein calculating Finvolves calculating an inverse Fourier transform of the algebraicexpression {tilde over (B)}(q)/{tilde over (h)}(q). wherein {tilde over(B)}′(q)is the Fourier transform of the image B′ of reduced noise, and{tilde over (h)}(q)is the Fourier transform of h.
 7. A method forprocessing a deconvoluted image B, the image B having been deconvolutedaccording to a deconvolution filter D, the method comprising reducingcorrelation between the image and the deconvolution filter, wherein anamount of correlation is calculated in a plurality of images P, and animage P′ is selected among the images P as an image having essentiallyminimal correlation with the deconvolution filter.
 8. The method ofclaim 7 wherein the amount of correlation C in an image P is calculatedaccording to the algebraic expression C=∫dq∥{tilde over(D)}(q)∥²·∥{tilde over (P)}(q)∥² wherein {tilde over (P)}(q) is theFourier transform of an image P.
 9. The method of claim 7 wherein {tildeover (P)}′(q), the Fourier transform of P′, is equal to {tilde over(P)}_(i+1)(q) for some integer i, where {tilde over (P)}_(i+1)(q) isobtained according to the algebraic expression {tilde over(P)}_(i+1)(q)={tilde over (P)}₁(q)(1+ε∥{tilde over (D)}(q)∥²)^(i), whereε is a small positive number.
 10. The method of claim 7 wherein {tildeover (P)}′(q) is obtained according to the algebraic expression e,otlP′(q)={tilde over (P)}₁(q)e^(−β∥{tilde over (D)}(q)∥) ² , where βis apredetermined constant.
 11. A method for processing a digital image B₁,the image B₁ being a convolution of an image F and a point spreadfunction h comprising the steps of: removing noise from the image B₁ soas to produce an image B′ of reduced noise; obtaining function {tildeover (P)}₁(q) according to the algebraic expression {tilde over(P)}₁(q)={tilde over (B)}′(q)/{tilde over (h)}(q); reducing correlationbetween {tilde over (P)}₁ and 1/{tilde over (h)} so as to product afunction {tilde over (P)}′ of reduced correlation; and obtaining arectified image F by inverse Fourier transform of {tilde over (P)}′(q).12. A method for obtaining a radius r of a point spread function hdescribing an out-of-focus distortion of a digital image B, the methodcomprising a step of calculating a gradient at a plurality of pixels inthe image B, in which a radius r(x) is calculated at each of theplurality of pixels based upon the gradient.
 13. The method according toclaim 12 wherein each of the plurality of pixels is located at an edgeof the image B.
 14. The method according to claim 12 wherein a radiusr(x) is inversely proportional to the gradient at x.
 15. The methodaccording to claim 13 wherein r is obtained as the r(x) having anessentially maximal frequency among the calculated radii r(x).
 16. Themethod according to claim 15 wherein a radius r(x) is calculatedaccording to the algebraic expression r(x)=2/πs(x) , wherein s(x) is theabsolute value of the gradient of B at x normalized by dividing by theheight of the edge at x.
 17. A method for processing a digital image B₁,the image B₁ being a convolution of an image F and a point spreadfunction h, comprising removing noise from the image B₁ so as to producean image B′ of reduced noise, and calculating F based upon B′, furthercomprising a step of producing the image B′ from an image B₀, where theimage B₀ was obtained using a digital camera that applies atransformation to a light level detected at a pixel, the transformationhaving an inverse, wherein B₁ is obtained from the image B₀ by applyingto the image B₀ the inverse transformation.
 18. A program storage devicereadable by machine, tangibly embodying a program of instructionsexecutable by the machine to perform method steps for processing adigital image B₁, the image B₁being a convolution of an image F and apoint spread function h, comprising the steps of: removing noise fromthe image B₁ so as to produce an image B′ of reduced noise; obtainingfunction {tilde over (P)}₁(q) according to the algebraic expression{tilde over (P)}₁(q)={tilde over (B)}′(q)/{tilde over (h)}(q); reducingcalculation between {tilde over (P)}₁ and {tilde over (1)}/h so as toproduct a function ′ of reduced correlation; and obtaining a rectifiedimage F by inverse Fourier transform of {tilde over (P)}′(q).